global tempoTeste

T = 0.1;
%% Maximos e minimos: nivel e sinal de controle
nivelMAX = 30;
nivelMIN = 0;
sinalContMAX = 3;%a1/Km*sqrt(2*g*nivelMAX);
sinalContMIN = -3;%a1/Km*sqrt(2*g*nivelMIN);

%% Range: nivel e sinal de controle
rangeNivel = nivelMAX - nivelMIN;
rangeCV = sinalContMAX - sinalContMIN;

%% Coletando dados
tempo = saidaMA(:,1);
CV = saidaMA(:,2);
tank1 = saidaMA(:,3);
tank2 = saidaMA(:,4);
ref = saidaMA(:,5);

%% FAIXA 1

t1 = tempo(tempoTeste/T:2*tempoTeste/T-1);
T2_1 = tank2(tempoTeste/T:2*tempoTeste/T-1);
CV_1 = CV(tempoTeste/T:2*tempoTeste/T-1);

deltaT2_1 = T2_1(end) - T2_1(1);
deltaCV_1 = CV_1(end) - CV_1(1);

K1 = deltaT2_1/deltaCV_1;

coef = ajpoly(t1',T2_1',size(T2_1,1),3);
x = t1';
y = coef(1) + coef(2)*x + coef(3)*x.^2 + coef(4)*x.^3;
fx = mdf(y,T);
fx2 = mdf(fx,T);

inflexao = min(find(fx2 < 0));
if inflexao < 2
   inflexao = 2; 
end
m = (y(inflexao+1)-y(inflexao-1))/(x(inflexao+1)-x(inflexao-1));
b = y(inflexao)-m*x(inflexao);

figure(1);clf(1)
plot(t1,T2_1,'b',x,y,'k',x(1:400),m*x(1:400)+b,'r',x(inflexao),y(inflexao),'.c')

theta1 = (T2_1(1)-b)/m - t1(1)
% tal1 = (T2_1(end)-b)/m - theta - t1(1);

aux = T2_1(1)+ (0.632*deltaT2_1);
index = min( find(T2_1 > aux) );
tal1 = t1(index)-t1(1)-theta1;

%% FAIXA 2

t2 = tempo(3*tempoTeste/T:4*tempoTeste/T-1);
T2_2 = tank2(3*tempoTeste/T:4*tempoTeste/T-1);
CV_2 = CV(3*tempoTeste/T:4*tempoTeste/T-1);

deltaT2_2 = T2_2(end) - T2_2(1);
deltaCV_2 = CV_2(end) - CV_2(1);
deltaT2_2_norm = deltaT2_2/rangeNivel;
deltaCV_2_norm = deltaCV_2/rangeCV;

K2 = deltaT2_2/deltaCV_2;

coef = ajpoly(t2',T2_2',size(T2_2,1),3);
x = t2';
y = coef(1) + coef(2)*x + coef(3)*x.^2 + coef(4)*x.^3;
fx = mdf(y,T);
fx2 = mdf(fx,T);

inflexao = min(find(fx2 > 0));
if inflexao < 2
   inflexao = 2; 
end
m = (y(inflexao+1)-y(inflexao-1))/(x(inflexao+1)-x(inflexao-1));
b = y(inflexao)-m*x(inflexao);

figure(2);clf(2)
plot(t2,T2_2,'b',x,y,'k',x(1:400),m*x(1:400)+b,'r',x(inflexao),y(inflexao),'.c')

theta2 = (T2_2(1)-b)/m - t2(1)

aux = T2_2(1)+ (0.632*deltaT2_2);
index = min( find(T2_2 < aux) );
tal2 = t2(index)-t2(1)-theta2;

%% FAIXA 3

t3 = tempo(5*tempoTeste/T:6*tempoTeste/T-1);
T2_3 = tank2(5*tempoTeste/T:6*tempoTeste/T-1);
CV_3 = CV(5*tempoTeste/T:6*tempoTeste/T-1);

deltaT2_3 = T2_3(end) - T2_3(1);
deltaCV_3 = CV_3(end) - CV_3(1);
deltaT2_3_norm = deltaT2_3/rangeNivel;
deltaCV_3_norm = deltaCV_3/rangeCV;

K3 = deltaT2_3/deltaCV_3;

coef = ajpoly(t3',T2_3',size(T2_3,1),3);
x = t3';
y = coef(1) + coef(2)*x + coef(3)*x.^2 + coef(4)*x.^3;
fx = mdf(y,T);
fx2 = mdf(fx,T);

inflexao = min(find(fx2 < 0));
if inflexao < 2
   inflexao = 2; 
end
m = (y(inflexao+1)-y(inflexao-1))/(x(inflexao+1)-x(inflexao-1));
b = y(inflexao)-m*x(inflexao);

figure(3);clf(3)
plot(t3,T2_3,'b',x,y,'k',x(1:400),m*x(1:400)+b,'r',x(inflexao),y(inflexao),'.c')

theta3 = (T2_3(1)-b)/m - t3(1)

aux = T2_3(1)+ (0.632*deltaT2_3);
index = min( find(T2_3 > aux) );
tal3 = t3(index)-t3(1)-theta3;

K = (K1+K2+K3)/3
theta = (theta1+theta2+theta3)/3
tal = (tal1+tal2+tal3)/3
fatorIncont = theta/tal

